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ABSTRACT 



A linear triangular finite element formulation was used 
to solve the steady state two dimensional conduction heat 
transfer equation. A FORTRAN IV computer program of the 
above formulation, employing double precision arithmetic 
and compact storage techniques, was applied to study heat 
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Numerical results of the heat transfer rate in the copper 
condenser were shown to have less than 2% variance to those 
obtained earlier. A rotating heat pipe designed to operate 
with a small value of outside heat transfer coefficient 
experiences no significant increase of the heat transfer rate 
with an increase in rotational speed, since the value of the 
outside thermal resistance dominates. 

Heat transfer rate continuously increases as the fin 
angle decreases. However, the increase is only slight when 
the fin angle is less than 11 degrees. 
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I. 



INTRODUCTION 



A. THE ROTATING HEAT PIPE 

The rotating heat pipe is a closed container designed 
to transfer a large amount of heat in rotating machinery. 

Its three main component parts are: a cylindrical evaporator, 

a truncated cone condenser, and a working fluid as shown 
in Figure 1. 

At rotation above the critical speed of a rotating heat 
pipe, the working fluid forms an annulus in the evaporator, 
and will be vaporized by heat addition in it. The vapor 
flows toward the condenser, as a result of a pressure differ- 
ence, transporting the latent heat of evaporation with it. 
External cooling of the condenser causes the vapor on the 
inner wall to condense and release its latent heat of evapor- 
ation. The centrifugal force due to the rotation has a 
component acting along the condenser wall that will act to 
drive the condensate back to the evaporator where the cycle 
is repeated. 

In a conventional heat pipe , the force driving the con- 
densate back to the evaporator is due to capillary action, 
which poses a limit to its operation. However, the driving 
force in the rotating heat pipe can be designed to operate 
in any orientation, including in a zero gravity field. 
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Figure 1 Schematic Drawing of a Rotating Heat Pipe 



B. UTILIZATION OF AN INTERNALLY FINNED CONDENSER 

The first theoretical investigation of the rotating 
heat pipe conducted at the Naval Postgraduate School, was 
performed by Ballback [1] in 1969. He studied the limi- 
tations of performance imposed on the rotating heat pipe 
due to various fluid dynamic mechanisms. Using existing 
theory and experimental correlations, he was able to esti- 
mate the sonic limit, boiling limit, entrainment limit, 
and the condensing limit of performance. 

1 . The Sonic Limit 

When increasing the heat flux in a rotating heat 
pipe, it is entirely conceivable to reach a limiting flow 
rate of the vapor brought on by a choked flow condition in 
the pipe. This condition imposes a limiting value on the 
amount of energy the vapor can transport, thus reducing the 
effectiveness of the heat pipe. Since the heat transfer 
method in the rotating heat pipe depends strongly on the 
latent energy of the working fluid, the rate of heat trans- 
fer in the heat pipe can be expressed as: 



Q 



t 



m 



v 




( 1 - 1 ) 



where m v = vapor mass flow rate in lbm/hr. From the con- 
tinuity equation, 



m 



' v 



= P 



v 



U 



( 1 - 2 ) 
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where 



U v = velocity of the vapor in ft/sec, and 

2 

A = cross sectional area for the vapor flow in ft . 
The heat transfer rate becomes 

Q. = p U A h* (1-3) 

*t v v fg 

and the vapor velocity is considered to be sonic, 

U v = c = v^kRT (1-4) 

where 

c = sonic velocity in ft/sec 

2 

gg = gravitational constant = 32.1739 f t-lbm/lbf-sec 
k = ratio of specific heats 
R = gas constant in ft-lbf/lbm R, and 
T = absolute temperature in °R. 

2 . Boiling Limit 

It has been postulated by Kutateladze [2] that the 
transition from nucleate to film boiling is totally a 
hydrodynamic process. With this postulation, he determined 
the theoretical formula for predicting the burnout heat 
flux: 
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K/ ^ A b h fg {a g(p f 



1/4 



(1-5) 



Q 



t 




where 



K 



constant value 



P v = density of the vapor in lbm/ft^ 

2 

A, = heat transfer area in the boiler in ft 



h^ = latent heat of vaporization in Btu/lbm 
a = surface tension in lbf/ft 

2 

g = acceleration of gravity in ft/hr 
= density of fluid in lbm/ft^ 



The experimental data obtained by Kutateladze suggested a 
value for K in the range of 0.13 to 0.19. 

3 . Entrainment Limit 

The flooding constraint in a wickless heat pipe was 
examined by Sakhuja [3], who assumed that: 



a) The vapor and liquid are exposed to each other 
in a counterflow movement. 

b) The counterflowing vapor tends to retard the 
falling film of liquid through shear stress. 

c) The liquid film remains stable and smooth, and 
the shear stress is usually small. 

d) The buoyancy forces resulting from the density 
difference between the vapor and liquid are the 
means for maintaining the counterflow. 
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e) The buoyancy forces are balanced by dissipative 
effects which are proportional to the momentum 
fluxes of vapor and liquid streams. 

His resulting correlation for flooding is 



Q 



t 



A x c2 h fq /9D(p f- p v )p v 

U + <p v /p f > 1/4 ] 2 



d-6) 



where 




heat transfer rate in Btu/hr 

2 

flow area in ft 

dimensionless constant, 0.725 for tube with 
sharp edged flange 

latent heat of vaporization in Btu/lbm 

2 

acceleration due to gravity in ft/hr 

inside diameter of heat pipe in ft 

3 

density of the fluid in lbm/ft 



P v = density of the vapor in lbm/ft 0 
4 . Condensing Limit 

Ballback [1] determined the condensation solution 
for a rotating heat pipe by modeling the condenser section 
of a rotating heat pipe as a rotating truncated cone, from 
which the following condensation limit was developed: 



, 2 k f 

"(I — 



“ 2 w w 3 * 



2 . 

sin (p 



} { [R + Lsin 4> ] - R } 



3 

4 



(1-7) 
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where 



Q 

k 

P 



t 

f 

f 



co 



fg 



w 



4 > 

R 



total heat transfer rate in Btu/hr 

thermal conductivity of the condensate film 
in Btu/hr-ft-F 

3 

density of fluid in lbm/ft 

angular velocity in 1/hr 

latent heat of vaporization in Btu/lbm 

saturation temperature in F 

inside wall temperature in F 

viscosity of fluid in lbm/ft-hr 

half cone angle in degrees 

minimum wall radius in ft 

length along the wall of the condenser in ft 



The condensing limit equation is a function of the geometry 
and speed of the rotating heat pipe, and the physical 
properties of the working fluid. 

Tantrakul [4] calculated these limitations on a 
heat pipe with specific physical characteristics as shown 
in Table 1, with the results shown in Figure 2. 

TABLE 1 

Specification of a Typical Rotating Heat Pipe 



Length 


14.000 


inches 


Minimum diameter 


2.000 


inches 


Wall thickness 


. 0.125 


inches 


Internal half angle 


1.000 


degree 


Rotating speed 


2700 


RPM 
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Heat Transfer Rate, Q Bt l u/hr 



2 

Saturation Pressure, P s lbf/in 




Figure 2 Operating Limits of a Typical , 

Water-filled, Rotating Heat Pipe 
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Obviously from the results in Figure 2, the condensing limit 
line is the predominant limitation for the amount of heat 
that can be transferred from the heat pipe. However, the 
other limitations may become important as the heat pipe 
geometry and operating conditions are varied. 

In order to augment the heat transfer capacity of 
the heat pipe, the current efforts are aimed at raising the 
condensing limit line which may be accomplished by: 

a) a high value of cone angle, to increase the 
driving force. 

b) some type of promoter of dropwise condensation 
to increase the value of the convective heat 
transfer coefficient, h. 

c) use of an internally finned condenser to increase 
the inner wall surface area and the value of 

h, since the presence of the fin will decrease 
the condensate film thickness. 

By machining fins in the inner wall of the condenser, 
a significant improvement of heat transfer will be realized. 

C. ANALYSIS OF THE INTERNALLY FINNED ROTATING HEAT PIPE 

A disadvantage of improving the performance of a rotating 
heat pipe with fins machined in the inner wall is that an 
increase in the wall thermal resistance results due to the 
increase in wall thickness brought on by the fins. By using 
a material with a high value of thermal conductivity, this 
effect can be minimized. 
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Figure 3 Internally Finned Condenser 

Geometry Showing Fins, Troughs 
and Lines of Symmetry 
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Figure 4 Internally Finned Condenser Geometry 
Considered in Analysis of Schafer [5] 
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Schafer [5], developed an analytical model for the case 
of fins with a triangular profile as shown in Figures 3 and 
4. In the analysis he assumed one dimensional heat conduc- 
tion through the wall. Corley [6] for this same case 
developed a two dimensional heat conduction model using a 

■ 

Finite Element Method, and also assumed a parabolic tem- 
perature distribution along the fin surface. His results 
indicated a significant improvement in heat transfer per- 
formance of about 75%, above that predicted by the one 
dimensional model of Schafer [5], However, Corley [6] 

cautiously noted a probable error of 50% existed in the 

\ 

calculation of the heat transfer at the fin -apex, and conse- 
quently mentioned that there may be a total heat transfer 
error of as high as 15%. Also, he noted that his analysis 
was not valid f£o a stainless steel condenser because of 
numerical difficulties. Tantrakul [4] modified Corley's 
computer program by increasing the number of finite elements 
in order to minimize the apex heat transfer error. His 
results with this modification converged with the results 
of Corley. 

D. THESIS OBJECTIVES 

The objectives of this thesis are: 

1. To generate a general computer program for two 
dimensional steady state Conduction Heat Transfer. 

2. To use this program to study heat transfer in a 
finned rotating heat pipe. 
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3. To conduct a parametric study for a proposed heat 
pipe design. 



II. FINITE ELEMENT SOLUTION 



A. REVIEW OF THE PREVIOUS ANALYSIS 

Schafer [5] studied the one dimensional model heat 
transfer solution and Corley [6] studied the two dimensional 
model for an internally finned rotating heat pipe. Both 
used the same assumptions and boundary conditions based 
upon the analysis of Ballback [1] , which are similar to 
those used in the Nusselt analysis of film condensation on 
a flat wall. The more important of those assumptions are: 

1. steady state operation 

2. film condensation, as opposed to dropwise condensation 

3. laminar condensate flow along both the fin and the 
trough 

4. static balance of forces within the condensate 

5. one dimensional conduction heat transfer through 
the film thickness (no convective heat transfer in 
the condensate film) 

6. no liquid - vapor interfacial shear forces 

7. no condensate subcooling 

8. zero heat flux boundary conditions on both sides of 
the wall section (symmetry conditions) , as shown 

in Figures 5 and 6 

9. saturation temperature at the fin apex 

10. zero film thickness at the fin apex 

11. negligible curvature of the wall 
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b. c. 



9 2 T 3 2 T n 

dx 7 ' dy 2 U 



a) T 1 ' T sat 



lj ) " k = h i ( T ' T sa t^ in ro 6 ion Cl] 

c ) ' k 377 = h 2 ( T " TJ in legion [ 3 ] 



<0 



3T 

3n 



= 0 in region [2] and [ 4 ] 



Figure 5 Differential Equation and Boundary 
Conditions Considered in Analysis 
of Corley [6] 
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Figure 6 Internally Finned Condenser 
Finite Element Geometry Con- 
sidered in Analysis of Corley [6] 
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The complete development of Schafer's finned condenser 
heat transfer is presented in his thesis. His results are 
used as a comparison in this thesis. 

However, since the same fluid flow and heat transfer 
theory applies to the two dimensional conduction model, 

Corley [6] utilized this model and further developed it 
using the Finite Element Method. He considered an internally 
finned geometry as shown in Figure 6, and used the fin 
condensate film equation derived by Schafer [5] 



6 (z) d6 (z) = 



k f (T g -T w (z) ) y f dz 



2 2 

p p oj (Rq + x sin 4>) h^ cos cos a 



(II-l) 



where 



6 = the fin condensate film thickness in ft' 

k f = thermal conductivity of working fluid in 

r Btu/hr-ft-F 

Pf = density of working fluid in lbm/ft^ 

Uf = viscosity of working fluid in lbm/ft-hr 

oj = angular velocity in 1/hr 

Rq = minimum radius in ft 

x = distance along apex of fin in ft 

Z = distance along the fin, from the apex in ft 

<p = half cone angle in degrees 

h^ = latent heat of vaporization in Btu/lbm 



He assumed a parabolic temperature distribution along the 
fin surface. 
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(II-2) 



V z) = r l z + s l z + fc l 



to calculate the convective heat transfer coefficient along 
the fin surface. He developed a solution of the steady 
state two dimensional conduction heat transfer equation: 



Vk • VT = 



( II- 3 ) 



with boundary conditions as shown in Figure 5 
a) saturation temperature at the apex. 



T = T 



sat 



b) convective boundary condition in region [1] 



-k 



3T 

3n 



h, (T - T ) , and in region [3] 

-L S3. "C 



-k 



3T 

3n 



h~ (T - t ) 
2 °° 



c) symmetry conditions, = 0, in region [2] and 
region [4 ] . 

He used the 8 node quadratic isoparametric element in the 
mesh shown in Figure 6, and a modification of the three 
dimensional isoparametric Finite Element Method computer 
program assembled by Lew [7] . 

Applying the boundary condition a) , T = T sa t at 2 = ^ 
gives t^ = T sat and equation (II- 2) becomes: 
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1 



or 



T ( z ) 
w 



r l Z 



S-^Z 



'sat ' 



\ 



(II-4 ) 



'sat 



- T w (z) = 



-r^z - s^z 



(II-5) 



rr - i 

Upon substituting equation (II-5) into equation (II-2) , 
and integrating the result from z = 0 to z , gives an 
equation for <5(z): 



6 (z) = [ 



4 k f (- ri \ \- )Vf 



1/4 



oj (R 0 + x sin <{)) h^ cos $ cos 



a 



✓ 



(II-6) 



f , \ 

Assuming one dimensional conduction heat transfer through 
the thin condensate film, the local convective heat transfer 
coefficient is: 



h (z) 



fin 



6 (z) 



= [- 



2 2 

a) (Rq + xsin cj>)h f ^cos 4>cos a 1/4 



4 M' r l 4 ~ 5 1T> 



-] 



(II-7) 

where k^ = fluid thermal conductivity in Btu/hr-ft-F which 
is applicable from z = 0 to z = Z* = Zq - 6*/cos a ; where 
6* is the trough condensate film thickness, as shown in 
Figure 5. Across the width of the trough, the convective 
heat transfer coefficient is: 
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trough 



(II-8) 



h(x) 



f 

5 * (x) 



The heat transfer rate for the entire condenser is then 
determined within the computer program algorithm by dividing 
its axial length into 100 increments of length Ax (see 
Figure 4) , and solving for the incremental heat transfer 
rate using an iterative procedure. The incremental rates 
are then summed to yield the total heat transfer rate. 

Figure 7 shows a flow chart description of the computer 

program used in Corley's analysis [6]. Constant values of 

r^ and s^ are determined at each increment using Lagrangian 

interpolation fit of the temperature at the nodes 1, 3, 

and 6 (see Figure 6) from the previous iteration. At the 

first increment, values of r^ and s n are calculated from an 

initial guess of the temperature at these nodes. The average 

nodal values of h, , hu and the average of h of the 

lavg 3avg ^ 

left hand boundary region associated with node 6 (see Figure 
6), are calculated using the position dependent value h(z)- in 



(equation (II-7)) in a Simpson's Rule approximation. In 
the trough, h( x ) troU gh is calculated from equation (II-3), 
using the 6* value calculated from the previous increment, 
and this value is used as the average convective heat trans- 
fer coefficient for nodes 9 and 12, and the right hand 
boundary of node 6. The average heat transfer coefficients 
calculated for the left- and right-hand boundary region of 
node 6 are weighted by their respective region length. These 
values are then added to yield the average convective 
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r 














Figure 7 Computer Program Flowchart of Two 
Dimensional Conduction Analysis of 
Corley [6] 
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The remaining boundary 



coefficient for node 6 , h 
conditions are now applied, and the incremental heat trans- 
fer rate per unit of condenser length Q i , and the nodal 
temperature values are generated from the Finite Element 
Method subroutine. If the resultant heat transfer rate 
converges to within a 0.05% change from the value of the 
previous iteration, i.e. 





< 0.0005 



(II-9) 



then the increment is deemed solved and is set equal to 
the final iteration, . However, if the convergence 
criterion is not met, newly calculated nodal temperature 
values are used to establish a new set of inner wall convec- 
tive heat transfer coef f icients , and the Finite Element 
solution is repeated until the heat transfer rate converges 
to within 0.0005 tolerance. The incremental mass flow rate 
in the trough is now calculated using the relation 



M 



tot i 




/! \ 

( 11 - 10 ) 



where this value is added to the mass flow from the previous 
increments. A modification of the trough mass flow rate 
equation from Schafer's thesis [5] is used to calculate the 
subsequent interval's trough condensate thickness <5*(x), 
with the polynomial root finder subroutine: 

U( \/ 
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M tot (x) 



2 2 2 
p w (R 0 + xsin 4>) 6* (x) sin (j> „ 

— [6 * (x) £ + 6 * (x) tan a} 



(II-ll) 



This incrementation is continued until the entire length of 
the condenser is transversed. Incremental heat transfer 
rates are summed, hence the total heat transfer rate: 

100 

Q, . = 2 N_p . J (Q.Ax) (II-12) 

tot fin L l 

i =1 

where N^ n is the number of fins along the inner wall 
circumference . 

To start the incrementation, two initial guesses are 
made. These are the initial value of 6* which is taken 
from the analysis by Sparrow and Gregg for condensation on 
a rotating disk [8] and the initial values of the tempera- 
tures at the nodes 1, 3 and 6. 

Corley [6] noted that an error of 50% due to calculation 
of heat transfer rate at the apex may cause a total error 
of 15%. To minimize this error he suggested an increase in 
the number of finite elements. 

Tantrakul [4] continued Corley's study. He increased 
the number of finite elements between the apex of the fin 
and the trough as shown in Figures 8 and 9 . To save com- 
puter execution time, he reduced the number of increments 
to 25 instead of 100 increments as in Corley's calculation. 
However, even with the reduced number of increments, the 
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t ** 0.0312S in. 
b - 0.025 in. 

R “ 0.762 in . 

I - 1 ° 



Figure 8 Internally Finned Condenser Finite 

Element Geometry (3 Finite Elements) 
Considered in Analysis of Tantrakul [4) 
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t - 0.03125 in. 
b « 0.025 in. 

R > ** 0.762 in. 

<*> - 1 0 . 

* ■ Co 



Figure 9 Internally Finned Condenser Finite 

Element Geometry (4 Finite Elements) 
Considered in Analysis of Tantrakul [4] 
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results of Tantrakul [4] converged favorably with Corley's 
results [6] . 

B . FORMULATION 

In this thesis the author used the same physical rela- 
tion and technique to calculate the heat transfer rate from 
the condenser. These are listed below: 

a) assume parabolic temperature distribution along the 
fin boundary 

b) utilize equation (II-6) to determine the condensate 
thickness , 6 (z) 

c) use equation (II-9), (11-10), and (11-11) to 

calculate the mass flow rate in the trough, and 
determine the subsequent condensate thickness. 

d) iterate to arrive at the desired solution with the 
convergence criterion 0.01, since the results are 
converged at this value (see Table 2) . 

For the Finite Element Method solution, the author used a 
linear triangular finite element model as shown in Figures 
10 through 12. 

Corley [6] solved the field equation for the fin with 
the following boundary conditions: 

a) temperature at the apex is equal to the saturation 
temperature 

b) convective boundary conditions in the inside and 
outside of the condenser 

c) symmetry condition on the sides. 
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I 




b = 0.025 in 
- 0.762 in 
0 = 1 ° 



Figure 10 Condenser Geometry Considered 
with 25 Linear Triangular 
Finite Elements 
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f 





/ 





t = 0.03125 in. 
b = 0.025 in. 

R 2 = 0.762 in. 
0 = 1 ° 



\ 

Figure 11 Condenser Geometry Considered 
with 40 Linear Triangular 
Finite Elements 
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/ y 



3 




C = 0.03125 in 
b = 0.025 in 

R 2 = °' 762 in < PV r 

0 - 1 ° 



Figure 12 Condenser Geometry Considered 
with 108 Linear Triangular 
Finite Elements 
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\ 



The first boundary condition indicates no heat transfer 
at the apex and a value of infinity for the convective heat 
transfer coefficient h. Since the value of h in reality 
never equals infinity, the author let the value of tempera- 
ture at the apex float, which is in agreement with the 
theoretical analysis for the extended surface in a one 
dimensional case. Only for a fin with a small angle or a 
very long length can the temperature at the apex become 
equal to the environment temperature (saturation temperature 
in this case) . 

Therefore the statement of the problem for the formula- 
tion of the Finite Element Method is shown in Figure 13: 

i? * = 0 (11-13) 

3x 3y 

with the boundary conditions : 

a) in region 1, -k -r- — = h n (T - T ) (11-14) 

on sat 

b) in region 3, -k (T - T^) (11-15) 

c) in region 2 and 4, = 0 (11-16) 



Dividing the domain into linear triangular finite elements, 
the differential equation within each element is: 



3 2 T (e) 3 2 T (e) 

3x 2 3y 2 



(11-17) 
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I 




i!l + i!l = 

8x 2 3y 2 



b.c. 



a ) 



b) 



c) 



- * Iff “ h i C T - T sat^ 1,1 rc S ion W 

' k 3n = h 2 (T " T J in rc e 101 ' W 

^ = 0 in region [2] and f-lj 



Figure 13 Differential Equation and Boundary 
Conditions Considered in Analysis 
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with the boundary conditions : 



a) 




3T (e) 
3n 




(11-18) 



for elements with sides along the convective boundary, and 



b) 



3T (e) 
3n 



0 



(11-19) 



for others. In the above equations. 



(e) 

(e) 



(e) 



T 

s 



(e) 



temperature within the element 

thermal conductivity of the element; assumed 
to be constant 

convective heat transfer coefficient along 
the boundary of the element, assumed to 
be constant 

surrounding temperature at the boundary of 
the element. 



Define the approximate value of the temperature distribution 
within each element as: 

3 

T (e) t T N (e) Cx,y) = l N i (x,y) T ± (e) 

1=1 ( 11 - 20 ) 

where 

= Two dimensional linear shape function, and 

( 0 ) 

= nodal temperature. 

Using equation (11-20) , equation (11-17) can be written in 
approximate form as 
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R = 



3 

I T i 

i=l 



( \ 3 2 N. 
(e) 1 

3x 2 



3 

I T i 

i=l 



/ . 3 2 N. 

(e) l 

3y 2 






where 



R 



= residual. 



Applying the Galerkin criterion: 



A 



(e) 



[R N. dA (e) = 0 

; 3 



( 11 - 22 ) 



and substituting equation (11-21) into equation (11-22) , 
gives : 



3 2 N. , . 

/I T i Ce) 3^ ^ ^ S 



(e) 



(e) 



, \ 3 2 N. , \ 

fjT. (e) — A. dA 6 = 0 

1 3y 2 3 



Applying Gauss' theorem to the factors on the right side 
yields the symmetric equation: 



, . 3N.3N. 3N.3N. , , 

/( ^ 2 + ^ArA)dA (e) 



A 



(e) 



3x 3x 3y 3y 



3T. 



(e) 



/N. — i dL (e) = 0 

1 1 3n 



(e) 



j 3 



(11-23) 
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M 



where 



" ( 0 ) 

A = area domain, and 

( 0 ) 

L = length domain. 

Substituting equation (11-18) and (11-19) into the second 
term of equation (11-23), gives: 



/ N 5 



3T . 

l 



(e) 



(e) 



3n 



dL 



(e) 



where 




<P (T (e) ) 



dL 



(e) 



(11-24) 



)Cf 21 ' 



K 

K 



/ 



<fi (T (e) ) 



(e) 

(e) 



(T 



(e) 



0 




for an element with its 
side along the convec- 
tive boundary 



others 



Therefore equation (11-23) can be written in the form: 

U- 
L 1 

i=l 



+ fix cj)(T (e) ) dL (e) = 0 (11-25) 



/[ 



3N. 3N • 
3 - 3 

3x 3x 



(e) 



3N . 3N . 
i 3 

3y 3y 



] dA 



(e) 
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The solution of 



/< 



(e) 



9N.9N. 3N.3N. , . 

i 1_ + 1 3 ) dA (e) 

9 x 3x 9y 3y ' 



(Pij ) 



(e) 



(11-26) 



where 



13 



b.b . 

1 JL 



c . c . 
1 9 



4A (e) 4A (e) 



(11-27) 



b. 

i 



Y i ' Y * 



(11-28) 



c i = X k - Xj (11-29) 

j \ = 1, 2, 3 
k 1 



as shown in Figure 14. Recalling equation (11-24), for 
elements with an edge along a convective boundary , 



(e) 



/Nj<J> (T (e) )dL (e) 



T. 

l 



(e) h 



(e) 



Te) 



(N.N. dL 
; l 3 



(e) 



(e) 






dL 



(e) 



(11-30) 



Applying equation (11-30) for elements with the nodal point 
1 and 2 located on the convective boundary, gives 
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3 




Figure 14 Linear Triangular Finite Element Geometry 
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(e) 



/ Nj <p (T (e) )dL (e) 



= T . 

l 



(e) h 



(e) 



(e) 



/ jP : 
(e) / 0 



0 | <p 1 ,p 2 ,0 > 

dL 



(e) 



- T 



(e) h 



(e) 



(e) 



/< P 9 } dL 



(e) 



(e) 



where 



(11-31) 



P 1 ' p 2 



= one dimensional linear shape function, 



Therefore, equation (11-23) can be written: 



(e) 



£t. (e) [ (p. . ) (e) + r 

L i *1] ,(e) 



(e) ( 0 



<p, p 9 0> 

L r Z , 



dL (e) ] 



h (e) T (e) 
k 1 * 7 s 



/ p 2 dL 



(e) 0 



(e) 



or 



(K) (e) {T} (e) = {F} (e) 



(II-31) 
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I, 3 



where 



) i 



«/ 



/ 

(K) 



(e) _ 



and 



r / 

2 2 

1 

+ - 


n (e) L (e) 
Q L 21 


. v4 

b l b 2 +c l c 2 , 


h <S, L 21 (e » 


1 

/ 

b-b^+c c. 
1 j 13 


4A (e) 


3k (e) 


4A (e) 


6k (e) 




b l b 2' h =l°2 


u ( e ) t (e) 
, h L 21 


, 2^ 2 h (e) (e) 

b 2 +c 2 ,h L 21 


/ 

b 2 b 3 +c 2 c 3 


4A (e) 


6k (e) 


4A (e) 


3k (e) 


4A (e) 


W c l c 3 


1-3 


/ 

b 2 b 3' te 2 c 3 




, 2/2 / 

b 3 4 c 3 y 


4A (e) 




4A^ e) 




4A {S) 










p 


(F> (e) 


h (e) T 


s (e) ^ < 


)[\ 


E3 




k <e) 1 


} o\ ' 


13 



3 



The element matrix equation is then assembled into a system 
matrix equation with the sequencing based upon element and 
system nodal point connectivity. 

The heat transfer rate within each element was calculated 



by 



rn (e) (e) 

Q (e) = h (e) L 12 (e > Ax(-i y— 2 T s (e) ) 



(11-32) 



where 



( 0 ) 

= temperature at node 1 



(e) 

T 2 = temperature at node 2 



<* 0 
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surrounding temperature 



s 




convective heat transfer coefficient 



L 



(e) 



length between nodal points 1 and 2 



21 



as shown in Figure 15 . 

C. THE COMPUTER PROGRAM 

The program consists of a main program and three 
subroutines ; 

a) the routine "COORD" used to define positions of 
system coordinate points 

b) the routine "FORMAF" used to formulate the Finite 
Element Method equations 

c) the routine "BANDEC" as an equation solver for a 
symmetric matrix which has been transformed into 
banded form. / 

The input data is entered with the units of inches for the 

2 

length, degrees for the angle, Btu/hr-ft -F for convective 
heat transfer coefficients, and in F for the temperatures. 

Input data is arranged in 6 classifications , and is 
shown in Table 3 . 
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Figure 15 • Heat Transfer Rate of an Element, Q 



TABLE 3 



Classification Of Input Data 

Classification Type of input data 

I system and element connectivity 

II condenser geometry 

III statement of the problem 

IV convergence criterion 

V specification to determine the system 

coordinate point 

VI initial guess of the temperature at the 

nodal points along the fin and the trough 

A detailed description of each of these classifications 
follows : 

Classification I 

READ (5, 100) NEL, NSNP, NBAN 



100 FORMAT (315) 






Explanation: NEL 


= number 


of elements 


NSNP 


= number 


of system nodal points 


NBAN 


= system 


band width 



READ (5, 200) (IEL, ( ICOR ( IEL , I ) , 1=1 , 3) , IEL=1 , NEL) 

200 FORMAT (.415) 

Explanation: IEL = the element number 

ICOR (IEL, I) = system nodal point (NP) corresponding 

to NP I of element IEL 



54 



Note: - numbering the element, first for all elements with 

a convective boundary, starting from the top to 
the bottom, then others 

- numbering the nodal point of an element is 
counterclockwise 

- the elements with convective- boundary have nodal 
point 1 and 2 located on the convective boundary 

Classification II 

READ (5, 300) CL, CANGL, RBASE , R2 , THICK, BFIN 
300 FORMAT ( 6G10.5) 

Explanation: CL = condenser length in inches 

CANGL = half cone angle in degrees 
RBASE = minimum radius of the condenser in 
inches 



R2 = RBASE / -- y - 1 

THICK = condenser wall thickness in inches 
BFIN = the height of the fin in inches 
READ (5, 400) NDIV, NEST, NEFB , NBOTI , NBOTF 



400 FORMAT ( 


515) 






Explanation : 


NDIV 


= number 
of the 


of increments along the length 
condenser 




NEST 


= number 
end of 


of the element on the right 
the trough 




NEFB 


= the element number with convective 
boundary located at the base of the 
fin 
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NBOTI 



the element number with convective 



boundary located at the right hand 
of the bottom side 

NBOTF = the element number with convective 
boundary located at the left hand 
of the bottom side 

Classification III 

READ (5, 500). (MRPM ( I) , 1=1 , 3 ) 

500 FORMAT ( 315 ) 

Explanation: MRPM = number of RPM ( co ) 

READ (5, 600) (TSS (I ) , 1=1 , 3) 

600 FORMAT ( 3G10.5 ) 

Explanation: TSS = saturation temperature of the 

working fluid ^ T sat^ 

READ (5, 700) TINF 
700 FORMAT (G10.5) 

Explanation: TINF = outside temperature (T^) 

READ (5, 800) (HINF ( I ) , 1=1 , 3 ) 

800 FORMAT ( 3G10.5 ) 

Explanation HINF = outside convective heat transfer 

coefficient (h . ) 

out 

Classification IV 
READ (5, 900) CRIT 
900 FORMAT ( G10.5) 

Explanation: CRIT = convergence criterion (0,01) 



56 



Classification V 



READ (5, 1000) (FANGL(I) , ( = 1,3) 

1000 FORMAT ( 3G10.5 ) 

Explanation: FANGL = fin half angle (a) 

READ (5, 1100) IFF 
1100 FORMAT ( 15 ) 

Explanation: • IFF = n-1, where n is the number of rows 

of the upper triangular fin section 
READ (5, 1200) (KFIN ( I ) , KFF (I ) , (=1 , IFF) 

1200 FORMAT ( 1615 ) 

Explanation KFIN = the number of system nodal points 

located on the symmetric boundary 
t—/ of triangular fin section, but does 
not include the system nodal points 
located at the base of the fin and 
the apex 

KFF = the number of system nodal points 
located along the fin convective 
boundary, but does not include the 
system nodal points located at 
the base of the fin and the apex 
READ (5, 1300) DOBF , DOTH, JTC, JLC , JINT, KT 
1300 FORMAT ( 2G10.5, 415 ) 

Explanation: DOBF = number of column within the fin 

DOTH = number of column within the trough 
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JTC 



= number of the system nodal point 
located at the junction of the 
symmetry boundary and the line of 
intersection between the fin and 
the condenser wall 
JLC = the number of the system nodal 
point located at the center of 
system coordinate 

JINT = the numerical difference between the 
two adjacent system nodal points 
vertically at the condenser wall 
section 

KT = the number of rows within the wall 
section 

Note : as an example for 25 finite elements is shown in 
Figure 16 
Classification VI 
READ (5, 1700) T(IE) 

READ (5, 1700) T(IG) 

1700 FORMAT CGI 0.5) 

Explanation: T = initial values of the temperature 

estimate at the nodal points 1 and 
2, for elements with convective 
boundary along the fin and the 
trough 
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Figure 16 Specification for Input Data to 
Determine the Coordinate of the 
System Nodal Points 
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D. TEST CASES 



Before the Finite Element Method program (routine 
"FORMAF" ) was applied to the condenser problem, a series of 
test cases were analysed. These cases allowed a comparison 
to the results obtained from previous analyses and reference 
9. These cases are listed below: 

1. A rectangular block in a uniformly convective 
environment on three sides, with a specified 
temperature on the non-convective boundary. These 
results are compared with reference 9 as shown in 
Figure 17. 

2. A rectangular block in a uniformly convective 
environment on the top and bottom side, insulated 
on the other sides. A comparison of these results 
with the exact values and Corley's results [6] is 
shown in Figure 18. 

3. A rectangular block with triangular fin in a uniformly 
convective environment along the fin and bottom side, 
insulated on the wall. A comparison of the results 
with Corley's results is shown in Figure 19. 

The results show good agreement with the results of reference 
9 (case 1) , and are comparable with the exact values (case 2) . 
Therefore, "FORMAF" can be used for further applications. 
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h = 120 Btu/hr-f t 2 - °F 




Note: numbers in parenthesis refer to the results of 

reference 9 



Figure 17 Statement and Results of Test Case 1 
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T = 100°F 

h = 100. Btu/hr-f t 2 -deg F 



Actual T = 24.1 °F 




T = 0 °F 

2 

h = 1000.0 Btu/hr-f t -deg F 



. 7586. Btu/hr actual 
Q = { 7552. Btu/hr 

1 (7527. Btu/hr) 



Note: Numbers in parentheses refer to results of 

Corley [6] . 



Figure 18. 



Statement and Results of Test Case 2 




c 

(ft) 


HEAT TRANSFER RATE (Btu/hr) 


ONE *) 

DIMENSIONAL 


TWO *) 

DIMENSIONAL 


TWO 

DIMENSIONAL 


0.25 


6526. 


6727. ' 


6826. 


0.05 


6561. 


6754 . 


6759. 



*) Results of Corley [6] 



Figure 19 Statement and Results of Test Case 3 



.0 ft 
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III. NUMERICAL RESULTS AND DISCUSSION 



The analysis concerns itself with the geometry of both 
a copper and stainless steel condenser as tabulated in 
Table 4, using water as the working -fluid. 



TABLE 4 

Specification of a Typical Internally 
Finned Rotating Heat Pipe 



Condenser length 


rz 


9.000 


inches 


Minimum radius 


= 


0.775 


inches 


Wall thickness 


- 


0.03125 


inches 


Height of fin 




0.025 


inches 


Cone half angle 


= 


1.0 


degree 


r J St ✓ a ' -cc 









Results are given in Tables 5 to 18 and Figures 20 to 28. 

It can be seen from the data that the solution of the Finite 
Element Method converges. This can be seen by examining 
data of temperature distribution within the fins of 25, 

40, and 108 linear triangular finite elements (Tables 7 to 
9), for the same location (nodes 10, 15, and 28 respectively). 

The heat transfer rate solution is generally affected 
in the same manner as the analysis of Schafer [5] when the 
outside convective heat transfer coefficient (h ) , rotational 
speed (oj) , fin angle (2a) , and thermal conductivity ( k wa j_j_) 
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TABLE 9 



Theoretical Heat Transfer Rate on a Stainless Stee^ 
Condenser at 1000 RPM, for h Qut = 1000. Btu/hr-ft 



Finned 

-deg.F 



Fin 

half 

angle 


No- 

of 

fins 


T „ 
sat 


Q(Btu/hr) 


(OF) 


25 EL 


40 EL 


108 EL 


10 


276 


100 


8948.2 


8944.1 


8893.9 






150 


23764.0 


23727.0 


23607.0 






200 


38557.0 


38492.0 


38286.0 


30 


84 


100 


8515.5 


8494.4 


8461.8 






150 


22511.0 


22445.0 


22354.0 






200 . 


36467.0 


36355.0 


36204.0 


50 


40 


100 


8186.3 


8164.6 


8141.6 






150 


21574.0 


21508.0 


21442.0 






200 


34916.0 


34805.0 


34693.0 
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TABLE 10 



Theoretical Heat Transfer Rate on a Stainless Steel Finned 

Condenser at 3000 RPM, for h ,_ = 1000. Btu/hr-ft 2 -deq.F 

out 3 



Fin 

half 

angle 


No. 

of 

fins 


T 

sat 

(°F) 


Q(Btu/hr) 


25 EL 


40 EL 


108 EL 


10 


276 


100. 


9096.4 


9100.4 


9058.1 






150. 


24193.0 


24200.0 


24098.0 






200. 


39277.0 


39284.0 


39123.0 


30 


■ 84 


100: 


8818.4 ' 


8810.3 


8786.0 






150. 


23389.0 


23376.0 


23301.0 






200. 


37934.0 


37900.0 


37781.0 


50. 


40 


100. 


8584.1 


8580.7 


8561.1 






150. 


22720.0 


22690.0 


22642.0 






200. 


36823.0 


36768.0 


36689.0 
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TABLE 11 



Theoretical Heat Transfer Rate on a Stainless Steel Finned 
Condenser at 1000 RPM, for h Qut = 5000. Btu/hr-f t 2 -deg . F 



Fin 

half 

angle 


No. 

of 

fins 


T 

sat 

(°F) 


Q(Btu/hr) 


25 EL 


40 EL 


108 EL 


10 


276 


100. 


26899.0 


26775.0 


26308. o' 






150. 


70520.0 




68673.0 






200. 


113950.0 


113300.0 


110780.0 


30 


84 


100. 


22742.0 


22540.0 


22286.0 






150. 


58879.0 


58288.0 


57582.0 






200, 


94734.0 


93748.0 


92581.0 


50 


40 


100. 


20426.0 


20253.0 


20086.0 






150. 


52480.0. 


51978.0 


51513.0 






200. 


84275.0 


83443.0 


82676.0 



TABLE 12 



■ i n pa f Transfer Rate on a Stainless Steel Finned 

Theoretical Heat Transfer = 500 0. Btu /hr-ft2-deg.F 

Condenser at 3000 KFM/ ro Q U t 



10 



30 



50 



276 



34 



40 



— r 

25 EL 1 


40 EL 


108 EL 


28705.0 


28701.0 


28263.0 


• 75738.0 


75734.0 


74492.0 


122610.0 


122620.0 


120570.0 


25604.0 


25534.0 


25300.0 


66822.0 


66734.0 


65951.0 


107830.0 


107800.0 


106260.0 


23673.0 


23555.0 


23410.0 


61488.0 


61127.0 


60726.0 


99069.0 


98464.0 


97803.0 
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TABLE 13 



Temperature Distribution Within the Fin at the 

10t“ Increment. 

(Copper Condenser with 25 Finite Elements 



NP 


T ( °F ) 


NP 


T ( °F) 


NP 


T ( °F ) 


1 


98.99 - 


15 


98.08 






2 


98.79 


16 


97 . 94 


* 




3 


98.77 


17 


98 .00 






4 


98.60 


18 


97.99 






5 


98.58 


19 


97.97 






6 


98.53 


20 


97.93 






7 


98.43 


21 


97.80 






8 


98.41 










9 


98.35 










10 


98.23 










11 


98.01 










12 


98.18 










13 


98 . 17 










14 


98 . 14 










a = 50 deg 




T = 70 °F 

oo 


^out 


= 1000 Btu/hr-ft 2 


-deg F RPM = 1000 


T sat : 


= 100°F 




Q 


= 9697 


. 2 Btu/hr 
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TABLE 14 



Temperature Distribution Within the Fin at the 

10 th increment 

(Copper Condenser with 40 Finite Elements) 



NP 


T(°F) 


NP 


T(°F) 


NP 


T ( °F) 


1 


99.01 


15 


98.22 


29 


97.90 


2 


98.84 


16 


98.04 


30 


97.82 


3 


98.83 


17 


97.99 


31 


97.79 


4 


98 . 69 


18 


98.17 






5 


98.68 


19 


98.16 






6 


98.64 


20 


98 . 14 






7 


98.55 


21 


98.11 






8 


98.54 


22 


98.06 






9 


98 . 50 


23 


97 .96 






10 


98.45 


24 


97.92 






11 


98.42 


25 


97.99 






12 


98.41 


26 


97.98 






13 


98 . 38 


27 


97.96 






14 


98.32 


•28 


97.94 






a = 5 


>0 deg 




T 


= 70 

CO 


°F 


h , 
out 


= 1000 Btu/hr-f t‘ 


'-deg F RPM = 1000 


T 


= 100 °F 




Q 


= 9687.2 Btu/hr 


sat 
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TABLE 15 



Temperature Distribution Within the Fin at the 

10th increment 

(Copper Condenser with 108 Finite Elements 



NP 


T(°F) 


NP 


T ( °F) 


NP 


T ( °F ) 


1 


99 . 04 


24 


98.39 


50 


97.90 


2 


98.91 


26 


98 .33 


51 


98.05 


3 


98 . 90 


28 


98 . 20 


52 


97.96 


4 


98.79 


29 


98 . 05 


58 


97.85 


6 


98 . 76 


30 


97.99 


60 


97.84 


7 


98.68 


31 


97.96 


61 


97.97 


8 


98.68 


32 


98.27 


62 


97.96 


10 


98.63 


34 


98.25 


63 


97.95 


11 


98.58 


36 


98.20 


64 


97.94 


13 


98.56 


38 


98.12 


65 


97.92 


15 


98.51 


39 


98.02 


66 


97.90 


16 


98.49 


40 


97.96 


67 


97.90 


18 


98.47 


41 


97.94 


68 


97.88 


19 


98.45 


42 


98.14 


69 


97.82 


21 


98.37 


44 


98.04 


70 


97.78 


22 


98.41 


48 


97 . 91 


71 


97 . 76 


a = 5 

h . : 
out 


0 deg 

= 1000 Btu/hr-ft 2 


T = 7 0 0 F 
00 

-deg F RPM = 1000 


T . : 

sat 


= 100 °F 




Q 


= 9678 


i . 9 Btu/hr 
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TABLE 16 



Temperature Distribution Within the Fin at the 

lO^-h Increment 

(Stainless Steel Condenser with 25 Finite Elements) 



NP 


T ( °F) 


NP 


T ( °F ) 


NP 


T ( °F) 


1 


99 . 79 


15 


95.17 






2 


99.08 


16 


94.44 






3 


99.20 


17 


93.94 






4 


99.19 


18 


93.91 






5 


98.24 


19 


93.83 






6 


98.37 


20 


93.69 






7 


97.23 


21 


93.11 






8 


97.21 










9 


97.10 










10 


96.66 










11 


95.48 










12 


95.50 










13 


95.47 










14 


95 . 36 










a = 50 


deg 




T 


= 70°F 


^out 


1000 But/hr-ft 2 - 


-deg F RPM = 1000 


Zsat 


100 °F 




Q : 


= 8186. 


, 3 Btu/hr 
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TABLE 17 

Temperature Distribution Within the Fin at the 

lO^h Increment 

(Stainless Steel Condenser with 40 Finite Elements) 



NP 


T (°F) 


NP 


T ( °F) 


NP 


T(°F) 


1 


99.82 


15 


96.70 


29 


93.62 


2 


99.31 


16 


95.73 


30 


93.26 


3 


99.39 


17 


95.50 


31 


93.12 


4 


98.67 


18 


95.50 






5 


98.72 


19 


95.48 






6 


98.84 


20 


95.41 






7 


97.97 


21 


95.30 






8 


97.99 


22 


95.12 








98 . 03 


23 


94.62 






10 


98.09 


24 


94.45 






11 


97.25 


25 


93.93 






12 


97.24 


' 26 


93.91 






13 


97.20 


27 


93 . 85 






14 


97 . 08 


28 


93 .75 






a = 50 


deg y 




T 


= 70 c 


'F 


h = 

out 


1000 Btu/hr-ft^- 


00 

•deg F y RPM = 1000 


T = 

sat 


100°F 




Q 


= 8164 


: . 6 Btu/hr 
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TABLE 18 



Temperature Distribution Within the Fin at the 

10 th Increment 

Stainless Steel Condenser with 108 Finite Elements) 



NP 


T ( °F ) 


NP 


T ( °F) 


NP 


T ( °F) 


1 


99.84 


24 


97.21 


50 


94.46 


2 


99.50 


26 


97.13 


51 


94 . 39 


3 


99.54 


28 


96 . 64 


52 


94.64 


4 


99.09 


29 


95.83 


58 


94.29 


6 


99.19 


30 


95.52 


60 


93.83 


7 


98.65 


31 


95.43 


61 


93.76 


8 


98.67 


32 


96.32 


62 


93.87 


10 


98.81 


34 


96.28 


63 


93.84 


11 


98.19 


36 


96.15 


64 


93.79 


13 


98.23 


38 


95 .82 


65 


93.73 


15 


98.36 


39 


95.31 


66 


93.64 


16 


97.70 


40 


95.03 


67 


93.55 


18 


97.72 


41 


94.95 


68 


93.31 


19 


97.73 


42 


95.45 


69 


93.31 


21 


97 . 73 


44 


95 . 41 


70 


93.13 


22 


97.22 


48 


95.04 


71 


93.07 


a = 50 


deg 




T 


= 7 0 °F 




h , = 
out 


1000 Btu/hr-f t 2 - 


CO 

-deg F RPM = 1000 


T = 

.. sat 


100 °F 




■ Q : 


= 8141. 


. 6 Btu/hr 
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Q (Btu/hr) 



200 



RPM = 1000 
a = 30° 

= copper 

- - - - = stainless steel 




< i i 

100 150 200 

Saturation Temperature (°F) 



2 

1. h = 5000 Btu/hr- ft - deg F 

out 

2 . h , = 1500 Btu/hr-ft 2 - deg F 

out 

2 

3. h = 1000 Btu/hr-ft -deg F 

out 

Figure 20 Heat Transfer Rate (Q) of Internally 

Finned Condenser at a Particular Value 
of h , vs. Saturation Temperature of 
out Working Fluid 
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2. a = 30° 

3. a = 50° 

Figure 21 Heat Transfer Rate (Q) of Internally Finned 
Condenser at a Particular Value of Fin Hair 
Angle vs. Saturation Temperature of Working Fluid 



80 



xie - ' 




3. h = 9000. Btu/hr-f t^-degF 



Figure 22 Heat Transfer Rate (Q) of Internally 
Finned Copper Condenser vs. RPM 



-degF 

-degF 
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3. ^ out = 9000* Btu/hr-f t 2 -degF 

Figure 23 Heat Transfer Rate (Q) of Internally Finned 

Stainless Steel Condenser vs. RPM 



-degF 

-degF 



32 



OO 3 



00a) 



001 



000 





No . of 






R. 

in 


Fins 


RPM 




1 


40 


1000 


T =100°F 


2 


200 


1000 


sat 


3 


40 


3000 


T ® =70°F 


4 


20 


3000 






■p'i i i i j i i i i j t r i' 

1C0C0 20000 30000 



i i r 
-10000 50000 



h (Btu/hr-ft -degF) 



Figure 24 



Thermal Resistance of 
Copper Condenser 



Internally 

vs . h . 
out 



Finned 
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003 



00c' — ; 



001 





No. of 






R. 

m 


Fins 


RPM 




i 


40 


1000 


T =100°F 


2 


200 


1000 


sat 


3 


40 


3000 


H 

8 

II 

o 

c 


4 


200 


3000 





000 




h OUt (Btu/hr ' ft 



Figure 25 Thermal Resistance of Internally Finned 
Stainless Steel Condenser vs. h 

out 
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a 



1. 


- = 0.01 
a 


4. 


1. 

a 


q f : 


heat transfer rate of 












internally finned con- 


2. 


- - 0.1 
a 


5 . 


£ = i. 

a 


5 


denser 


3. 


- = 0.5 


6. 


c ~ 

— « 2 . 


% : 


heat transfer rate of 




a 




a 




smooth condenser 




Figure 27 


Compar 


ison of 


Heat Transfer Rate (Q w /Q q ) 






b 






2 






vs . — 
a 


at h , 
out 


= 1000. 


Btu/hr-ft -degF 
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1 

2 . 

3 . 



- » 0.01 

a 



- “ 0.1 

a 



- =■ 0.5 

a 



4 , 

5 , 

6 , 



= l. 



- *> 1.5 
a 



- - 2 , 



Qp : heat transfer rate 0 f 
internally finned con- 
denser. 

Qg : heat transfer rate of 
smooth condenser. 



Figure 28 Comparison of Heat Transfer Rate (Q„/Q c ) 

b 2 F S 

vs. —at h = 5000. Btu/hr-ft -degF 
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were varied. In Figure 20, the influence of the outside 
convective heat transfer coefficient, h , is shown both 
for copper and stainless steel. Figure 21 illustrates how 
the rate of heat transfer was affected when the fin half 
angle was varied. In addition, the heat transfer rate was 
found not to vary significantly with rotational speed when 
the heat pipe operated with low values of outside convective 
heat transfer coefficient, h . However, for high values 
of the effect of rotational speed cannot be ignored 

(see Figures 22 and 23) . 

The small variance of the heat transfer rate with rota- 
tional speed is due to the dominancy of the outside thermal 
resistance, R ou1 _/ such that changes in the condensate film 
thickness, 6, will not significantly affect the total thermal 
resistance as shown in Figures 24 and 25. The heat transfer 
rate will have a large variance when the heat pipe is operated 

with a high value of h 

’ out 

From the results of the parametric study, as shown in 
Figures 27 and 28, the heat transfer rate continuously in- 
creases as the fin half angle decreases (b/a increases). 

Since the total number of the fins in the condenser is 
larger with a small half angle fin, this results in an in- 
crease in the heat transfer rate. However, the increase 
is only slight when the fin half angle less than 11 degrees 
(b/a greater than 5) . 

The study shows the enhancement in heat transfer is 
larger when operating with a high value of or when 
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Rout is small. Therefore the design operating with a 
high value of h using both internal and external fins 

is recommended. 
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IV. CONCLUSIONS 



1. The Finite Element Method works and converges. 

2. The designer has a choice of: 

- material 

- number of fins and fin angle 

- RPM 

- h 

out 

3. The heat transfer rate increases continuously 
with decreasing fin angle; however, for half 
angle less than 11 degrees, the increase is 
only slight. 
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V. RECOMMENDATIONS 



1. Manufacture internally finned heat pipe and test. 

2. Manufacture internally and externally finned heat 
pipe and test. 
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APPENDIX 



COMPUTER PROGRAM LISTING 



*******X*XXXX**r*X*4 X*******************;*-**.**** 

X 

X 

TWO DIMENSIONAL heat conduction finite element analysis X 

OF ROTATING HEAT PIPE , USING TRIANGULAR ELEMENT MODEL * 
COMPILED FV MAJOR IGNATIUS . S . PURNOMO IN JUNE ,1973 X 

X 

t 

XXXXXtt X'VXXXXXXXXTXXXXXXXXXXXXXxXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
IMPLICIT PEALX3(A-H,0-Z) 

DIMENSION Z < 200 ) , EPS ( SOU ) , HZ < 200 ) , XC0FC5 ) . COF < 5 ) . RQOTr < 4 ) , ROOT I ( 
t . T ( 200 ) . GE C 200 ) , DMDOT HOC ) , Ur ( 200 ) , CF ( 20O ) , RHOF ( 200 ) , DEL ( 203 ) , CU 
tO0 ) • AMTOT r 200 ) , R ( 200 ) • QINC ( 200 > , TB 1 200 ) , TT ( 200 ) , TEN C 200 ) , TE ( 200 ) 
tIBf 200 ) • NO ( 200 ) 

DIMENS I Of I T9S <. 3 ' , FANGl C 3 ) . HINF (3 0 MRPMC 3 ) 

COMMON/ APOL/DOBF , DOTH , KF IN( SO ) , KFF ( SO ) , IFF . JTC , JLC , J INT , KT 
COMMON/MAr 0/A ( 200 - 50 ) , F ( 200 ,U,H( 200 ) , TS(2O0 ) . TSAT , CK , NEL , NSNP . M 

*N. I COR ( 200 ■ 3 ) 

COMMON/PCRD/X ( 200 ) , Y< 200 ) , EZEPO, BUIN , THICK , TALFA , APS 
HDEN ( A1 , B1 , ZZ)-(-l . 0D0*(AUZZ**3/3.0D0+BUZZ*X2/2 0DO ) ) 

ELEMENT CONNECTIUITIES 

READ(S, 100 ) NEL . NSNP , MS AN 
O FORMAT (3 IS) 

UR I TEC 6, ISO) NEL, NSNP > MEAN 

') FORMAT ( /2X , 'NO . OF . ELEMENTS" ' . IS, 10X . ' NO . OF SYSTEM N.P -'.IS.10X. 



XO OF EANDED= ' , IS) 



READ 5 - 200 ) ( IEL, ( ICORf I EL, I ), 1-1,3). IEL = 1 , NEL ) 

3 FORMuT ( ^IS ) 

URITE C 6 - 250 ) 

5 FORMAT ■' .-'2X • ' ELEMENT ‘ , 1CX, 'NP1' , 14X, 'UPS’ , 15X, 'NP3' ) 
URITE '6.251 KIEL, cICORdEL, I ), 1-1.3). IEL-l.NEL) 

5 FORMAT ( IS. 12X. IS. 12X. IS. 12X. IS) 



THE CONDENSER GEOMETRY 

CL. IS THE CONDENSER LENGTH 

CANGL- IS THE CONE HALF ANGLE 

RBASE IS THE CONDENSER RADIUS AT THE BASE 

THICK IS THE WALL THICKNESS OF THE CONDENSER 

READ f S • 300 ) CL . CANGL . RBASE , R2 , THICK . EFIN 
30 FORMAT < 6G10 5) 

CL=CL''12 ODO 
R2=R2/'12 0 DO " 

RBASE = RBASE/ 12 ODO 
B"IN=BFIN/12 ODO 

' NDIU IS THE NUMBER OF INCREMENT 
NEST IS THE ELEMENT NUMEER AT THE END OF THE TROUGH 
NEFB IS THE ELEMENT NUMBER AT THE BASE OF THE FIN 
NBOTI IS THE FIRST ELEMENT AT THE BOTTOM SIDE 
NBOTF IS THE LAST ELEMENT AT THE BOTTOM SIDE 

READ(S. 400 ) NDIO. NEST, NEFB. NBOTI, NBOTF 
•0 FORMAT (5 IS) 

DIU = DFLOAT ( NDIU ) 

PI *3 1 41592E3353979D0 
PHI =2 0D0*CANGL*PI/3S0 OD0 
SPHI»DSIN(PHI ) 

CPHI“DCOS(PHI ) 

TFH I * DT AN (PHI ) 
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DELX-CL^Dlu 

CEASE* 2 . ODO*PI XRBASE 
REXIT-PBA5E+CLXTPHI 
CEXIT-2 ,0D0*PI*PEXIT 
THICK-THICK/12.0DQ 

URITE f 6, 450)CL,CANGL .R2. THICK .RErtSE, EUIN 
30 FORMAT*.' ' , /V5X, 'CONDENSERS LENGTH- ' ,E12 . 5. 30X, 'HALF CONE ANGLE- 

#,E12 5, s.SX. 'R2-' ,E12.S,45X, 'THICK- ' ,E12 . S, x.SX. 'R-BASE- ' , E12 . 5 . - 

XX, 'B-FIN- ' . E12 5) 

DATA FOR RUNNING 

READ < 5 • 500 ) (MRPMCI ), 1-1,3 ) 

•50 FORMAT < 315 ) - 

READ (5, 600 ) (TSSCI ). 1-1,3) 
i JO FORMAT C.3G10 5) 

PEpOC 5 . 700 ) TINF 
00 FORMAT CG10 5) 

READ (5 , 800 > ( HINF (I). 1-1,3) 

30 FORMAT (3G1 0.5) 

THE CON'JERGENCE CRITERIAN 

READ ( 5 , 900 ) CRIT 
30 FORMAT (G10 9) 

INTERNALLY FIN GEOMETRY 

READCS, 1000 ) (FANGL(I), 1-1,3) . 

30 FORMAT ( 3G10 . 5 ) 

READCS. 1 1O0 ) IFF 
30 FORMAT CIS) 

w "G 
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/ 











RErtC'5/ 1200) (KFIIN I ),IIFF(I ) . I * 1 ■ IFF) 

20 FORMAT ( 16 IE ) 

READ*; 5, 1300 )D0BF , DOTH . JTC , JLC , JINT . KT 
30 FORMAT » EG 10 .5,415) 

MHB-HEFBX2 

' NTM = NB0TI-KN'60TF-NB0TI V2 
MEF =»f fDOTF+1 

tJRITE ( 6 , 1350 ) ICORCNBOTI , 2 ) , ICORCNEFB, 1 ) - ICORCNTM, 2 ) - ICORCNEST. 
*ICOR(NBOTF.l ) 

33 FORPIAT ( ///5X, 'TIB*' , 15, 10X, 'TT = ' , 15, /.SX, 'TBfT= ^ , 15, 10X, 'TE = ' , IS. a 

X6X. 'TB*' • 15) 

DO 1-400 KIh-1.3 

ALFA-FANGLUCIA >*2 0D0XPI/360 ODO ' 

SALFA-DSIMALFA) 

CALFA-DCOS(ALFA) 

TALFA*DTANi ALFA) 

EZER0 3 2 ODOXBUINXTALFA 
EPSO=EZERO 

NFIM-CBASE^CEZERO+EPSO ) 

EPSEX" ( CEXIT- ( NFINXEZERO ) )/NFIN 
BETA* < EP5EX-EPS0 )/DIV 
ZZERO-BUIN/CALFh 
ZA 2 0 . CD0 

BOUNDARY CONDITIONS AND TEMPERATURE'S ESTIMATE 
ALONG THE FIN BOUNDARY 

DO 1450 NT I Nr =NBCT I , N30TF 
.3 TSCNTINF )-TINF. 

DO 1500 NNT-NSF.NEL 
TS ( NMT ) »0 0D0 x 
;3 H ( NNT ) • 0 0D0 



> 



i 
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DO 1600 IGT-i,MEST 
IE-IC0R(IC»T,2> 
t'O REhD‘5, 1700) T(IE) 

1G“ ICOP. ( NEST - 1 i 
PErtt'tS. 1700) TCIG) 

"iO FORMAT (G 10 5) 

DO mOO KIB-1.3 

nRPH=nppruK;ip ) 

OMEGA -NRPMI'2 0DOTPIX60 ODO 
DO 1400 KIC-1,3 
DO 1200 KL-MBOTI.NBOTF 
S*C HC OL )=*HINF (KIC ) 
hTfm^hinf (KCIC ) 

DO 1400 r.ID-1,3 
T5fiT-TS5(KID) 

DO 1900 MS AT »l. NEST 
90 TS'NSAT >«TShT 

TSOLID= (T5A7+TIMF )/2 . ODO 

DEUD-0 OOOOG7F2DO 

QT-J ODO 

OTOT =0 ODO 

DMTOT-0 ODO 

NK-NDIU+1 

DO 2000 NI-l.NK 

R( NT )»R2+NIXDELX*SPHI 

EPSCNI )-EPSO+NI*BETA 

APS»EPS(NI ) 

NODAL POINTS COORDINATE 

CALL COORD 
2(1 ) *ZA 

DO 2100 IZEL-l.NEFB 
NA-ICORC IZEL, 1 ) 
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(SB • ICOR >. IZEL , 2 ) 

XE-X«HA , '~X( NO ) 

VE«Y*NA )-Y<NB ) 

ELZ-DSGPT < X£T + 2+VE;U2 • 

!1)0 ZiIZCL+1 ) = Z< IZED+ELZ 

XZB«X».ICOR(NHB • 1 ) )-X(ICOR< 1 .2) ) 
S'ZB-Y < ICOR* NHB, 1 * )-Y< IOORC 1.2)) 
ZE=BSGPT < XZB*4:2+VZB**2 ) 
ZC-ZZEPO 
. IN-1 



FARAEOLIC TENPERhTURE DISTRIBUTION ALONG the fin 
BOUNDARY. USING LAGRANGE INTERPOLATION 

> TP1-TC ICORCl .2 ) ) 

TP2-T* ICOR( NHB . 1 > ) 

TP3-T c ICOR (HEFB , 1 ) ) 

AP1=TP1/(ZB*ZC ) 

AP2 = TP2/(ZB4*(ZB-ZC ) ) 

AP3»TP3/(ZCX(ZC-ZB ) ) 

EP1 =-( ZB + ZC UAP1 
BP2--ZCXRPH 
BP3= — ZBX AP3 
A1-AP1+AP2+AP3 
B1-BP1+BP2+BP3 
TC-0 0DO 

DO 2200 NY- 1. NEST 
3)0 TC=TC+T( ICCRCNY, 2) ) 

AY-DFLOATCNY*! ) 

TF «(TC-t-TClCCR( NY . 1 ) )+AYXTS( NY ) U C2 . ODOXAY ) 

SOLID-FLUID PROPERTIES 
HFG-1097 . 2D0-0 . 60187SDQ*TS(1 ) 
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RHOF'NI 1-62 . 774D0-0 00255638D0*TF-0 000O53572D0*TF*#2 

CF '.HI ) = 0 . 3034D0+0 OQ073S927DQXTF-0 OO00O147321D0*TF**2 

UFCfll )«0. O01397DO-0. 00001 4669P0*TF+O.OO000006312S3D0*TFX*2-0.000~ 

) 

*O0Q0097G569DO*TF**3 
UF (HI )-3600*UF<HI ) 

CU cm ) = 231 . 777ED0-0 02222D0XTS0LID 
CK-CUKMI ) 

CONST a RHOF (MI ) T YSTOnEGAXXEXHF GXCPH 1 1 C A LF AXR ( N I ) 

AVERAGE ELEMENT CONVECTIVE COEFFICIENT ALONG THE 
FIN BOUNDARY 

ZSTAR a ZZERO-DEL C N I )/CALFA 
AZZ-DELCNI )/SALFA 
ZZ = 7, STAR 

AZS=DABS( 4*CF ( NI )*UF ( NI )XHDEN ( A1 , B1 . ZZ )/CONST )X*0 25D0 

HAC=0 ODO 

DO 230O IEL-l.NEFB 

AZ*Z( IEL ) 

BZ-ZCIEL+1 ) 

IFCZSTAR LE.BZ) GO TO 2400 
GO TO 2500 

'•50 IFCHAC NE 0 ODO) GO TO 2600 
BZ=Z9TAR 

5)0 IF (IEL NE 1 ) GO TO 2700 
AK=(8Z-AZ)/S.0D0 
ZZ-AK 

GO TO 2800 
!‘)0 A<* ( BZ-AZ )/4 . 000 
ZZ-AZ 

use zel=4*ak 

DO 2900 NH*1 .5 

HZCNH )■ DABS (OF (NI )#*3*C0NST/(4*UF(NI )*HDEM( A1 , Bl , ZZ ) ) )XX0 25D0 
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3$)0 ZZ-ZZ+AK 

CONH-Ata CHZ( 1 >+4*HZ<2)+e*HZ<3 )+4*HZ(4 )+HZ(S ) )/(3*Z£L ) 
IFCZSTAR EQ.BZ) GO TO 3000 
H( IEL )-CONH 
GO TO 2300 
300 AZ-ZSTAR 

HAZ-CONH*' AZ~Z< IZL ) ) 

DELfi-AZS 
300 BZ-Z<IEL+1> 

DELB= ( BZ-ZSTAR )f AZZ/ (ZZERO-ZSTAR ) 

DhLZ a f. BELA+DELB V£ ODO 
HAC- ( BZ-AZ )£CF ( NI )/DELZ 
H( IEL )• (HAZ+HAC )/(BZ-Z( IEL ) ) 

GO TO 2300 
300 AZ-Zf IEL l 
BELA* DELB 
HrtZ =>0.000 
GO TO 3100 
300 CONTINUE 

NETI-NEFB+1 
DO 3200 I EL “NET I . NEST 
300 H(IEL)“CF(NI VDELCNI ) 

ENTRY INTO THE FINITE ELEMENT SOLUTION 

CALL FORMAF 

CALL BANDEC ( NSNP • NBAN • 1 ) 

THE TEMPERATURE DISTRIBUTION 

DO 3300 NT- 1 -NSNP 
100 T(NT ) -F (NT, 1 ) 

TIB ( NI )-T< ICORCNBOTI .2) ) 

TT(NI ) *T ( I COR ( NEFB . 1 ) ) 
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TBPKNI >«T< ICOPCNTN. 2 ) ) 
TE ( ft I ' • T f ICOR' NEST, 1 ) ) 

te<ni ,<*t< icorcnbotf, i ) ) 

TTS-0 0D0 
DO 3-100 NS 3 1 , NSNP 
}-)0 TTS-TTS+T ( NS ) 

PN-DFLOAT '.NS > 
TSOLID-TTS/PN 



Q AT THE BOTTOM SIDE 
QE I 3 0 ODO 

DO 3S00 IBEL-NBOTI.HBOTF 
NKA-ICORt IFEL, 1 > 

NKB-ICQR(IBEL.2> 

XB 3 X<NKA J-XCMKB ) 

YB-YCNKA J-YtNKB ) 

ELB 3 DSQRT(XB**2+Y3**2 ) 

3*50 QBI-QEI + (T(MKA)+T(NKB)-2*T5<IBEL))*ELB*H(IBEL)X2 0D0 
QB^NI )-QBI*DELX 



ITERATION UNTIL CONVERGENCE CRITERIA IS NET 

IF(IP1 EO 1 ; GO TO 2600 
QJ-QBI 
GO TO 3700 
3130 QI 3 QBI 
IH-2 
GO TO 2 

T30 A0>DA3S(QJ-QI )/0J 

IF (AO LE.CRIT) GO TO 3800 

QI-GJ 

GO TO 2 

130 DMDOTCNI )«2*Q3IXDELX/HFG 
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DP1T0T -Dl'iTOT +DMDOT • HI ) 

Cl«RHOF(HI )**2*0P1EGA**2*R(MI )/(3*UF(Nl ) ) 

XCOFC 1 >--DMTOT 
XCOFi‘2 -'>0 OD0 
XCOF ( 3 ) *0 0D0 
XCOFM >-CiTEPS f Ml ) 

XCOF ( 5 ) - CITThLFA 
f1 a 4 

CALL DF'OLPT ( XCOF . COF , H , ROOTR , ROOT I . IER ) 

IF ( R007R ( 1 ) . LT 0 C05D0 AMD RCCTR(l) GT 0 0D0 ) GO TO 2900 

IF C ROOTR r 2; . LT 0 OOSD0 AMD R00TRC2) GT 0 . 0D0 ) GO TO 4O00 

IF ( ROOTR < 3 ) LT 0 COSDO \ND RCCTRC3) GT.0 CD0 ) GO TO 4100 

IF ( ROOTR r A V LT. O OOSDO AMD ROOTR ( 4 ) GT . 0 . 0D0 ) GO TO 4200 

URITEfb 4300) 

:>0 FORMAT ( ''/10X , ' CRASH , CRASH . CRASH ' ) 

URITECb, 4350 ) ( ROOTR ( I ) , I* 1 , 4 ) 

C!0 FORMAT (//5X. 4 C El 2 7 . 3X ) ) 

GO TO 61O0 

THE CONDENSATE THICKNESS 

00 DEL C N I + 1 ) -ROOTR ( 1 ) 

GO TO 4400 

•00 DELCMI+l ) -ROOTR C 2 ) 

GO TO 4400 

1) 0 DELCNI + l J-ROOTRO ) 

GO TO 4400 

2) 0 DELCNI + l )=R00TR(4 ) 

00 QEL s O 0D0 

IFCNI ME 10) GO TO 4S00 
(JRITE ( S . 4600 ) 

00 FORMAT C ' 1 ' - //SX, ' NR ' , 6X , ' X ' , 12X , ' Y 7 , 12X, 'T' ) 

DO 4700 NP-l.NSNP 

’,0 URITEC6. 4300 ) NP , X ( NP ) . VCNP ) , T (NP ) 
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IS'O FORMAT (/ax, I3,3X.3<F10.6.3X)) 

URITE)' 6. -4900) 

FORMAT (K/2X, 'EL' ,3X. 'H MIX, 'EL-LENGTH' , 15X, 'Q-EL' ) 

DO 5000 KKL»1,NB0TF 
HKX= ICORC KKL, 1 ) 

nky»icor(kkl.2) 

XP=X« NKX )-X<NKY ) 

YP-Y<NKX )-Y(NKY ) 

£XY*DSQRT < XP**2+YP* *2) 

QEP • Drtfr'5 r ( T ( NKX )+T ( NKY )-2*TS ( KKL ) )*EXY*H < KKL ) /2 . ODO ) 
GEP a QEP*DELX 

JOC URITECG, 5100 ) KKL , H < KKL ' , EX Y . OEP 
3)0 FORMAT < V2X. 12, 3X,E12 5,3X,E12 S.10X.E125) 

URITE r <5 • 5200 ) CRIT 

,'c)0 FORMAT < x2X. 'CONUERGENCE CRITERIAN- ' , E15 . 8 ) 

Q FROM THE TOP SIDE 

♦5)0 DO 5300 IQEL-l.NEST 
KA-ICOR(IGEL.l > 

KB* I COR C I GEL- 2 ) 

XQEL"X(KA )-X( KB ) 

YOEL"Y(KB )— Y ( KA ) 

ELM" DSGRT ( XQEL**2+YGEL**2 ) 

QEL-QEL + ( 2*TS ( IGEL )-T ( KA )-T ( KB ) )*ELM*H< IOEL )/2 . ODO 
>1)0 CONTINUE 

QINCCNI )»GEL*DELX 
AMTOTCNI ) "DNTOT 
GET »QELxDELX*NFlN*2 
QT »QT +GET 

GA-GBri'DELXXNFIN*2 
QTOT-QTCT+GA 
50 CONTINUE 

URITE C 6 . 5400 ) HFG . NFIN - H ( NBOTF ) . T5AT . NRPM , QTOT . QT , FANGL ( KI A ) 
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> 

+0 FORMATS ' ' ,/,',5X, 'HFG-' ,E12 5, /, 5X, ''NUMBER OF FIN- ' , 15. / , SX , 'H~0~ 

. E12 S.x, SX. 'TS AT-' . E12.S, /.SX. 'RPN»' , IS, /,5X. 'G-BOT- ' ,E12.5, 

> | 

*X, 'Q-TOP 3 ' . E12 . 5 , / , 5X , 'HALF-ANGLE 3 7 , F3 3 ) 

WRITE re, 5500) 

5*0 FORNAT ( 'O' ,6X, 'J' ,4X, 'FILN THICKNESS' ,6X, 'Q-INCPEM' . 6X, '(1ASS-TCT~ 

17X. 'TIB ' ,3X- 'TT' , 10X, 'TE' ,2X. 'TB' ) 

DO 5600 NR-l.NDIO 

€i0 WRITE < S • 555 ) NR , DEL ( NR ) . OB ( NR ) , AMTOT ( NR ) , TIB C NR ) , TT ( NR ) , TE ( NR ) . T~ 
xriR) 

V FORNAT (' ' .4X, I4,4X,F12 19, 4X, F10 . 4 , 6X, F9 5 , 6X , F5 . 1 , 6X. F5 . 1 , 6X, ~ 

1 . 1,6X,F5 1 ) 

WRITE (6. 5309 ) 

<£?0 FORNAT ( '0' , 6X, ' J' ,6X, 'K-WALL' , 4X, 'K-FILN' ,3X, 'DENSITY' , 4X, 'UISC-- > 

*ILN' . 6X . 'EPSILON' ,5X, 'RADIUS' ,5X, 'TEN' , 5X, 'Q-BOT' ) 

DO 5900 NG 3 1 , NDIU . 2 

;<)0 UR ITE ( 6 . 777 ) NG . CU < NG ) , OF ( NG ) , RHOF ( NG ) , UF ( NG ) , EPS ( NG ) , R ( NG ) . TBfl < - 

I 

* ) . QINCCNG ) 

•' FORNAT (' ' ,4X, I4.4X.F7.3,4X,F6.4,4X,F6.3,4X,F9.?,4X,F9.7,4X,F7 ■*» 

15X, F5 1 , IX, 1P1D12 3) 

00 CONTINUE 
00 STOP 
END 

SUBROUTINE COORD 
INPLICIT REAL*8(A-H.0-Z) 

COmON/PCRD/X ( 200 ) . Y C 200 ) - EZERO , BUIN . THICK , TALFA , APS 
COnnCN/APOL/DOBF , DOTH , KFINC50 ) , KFF ( 50 ) , IFF . JTC - JLC . JIMT . KT 
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DELH-BUIN'DOBF 
X( 1 ) *0 OD0 
Yi, 1 ) =TH1CK+BUIN 
M-l 

DO 3£1 1 = 1 - IFF 
ICA-KFINC I ) 

ICB-KFF« I ) 

CBA-DFLOrtT(ICB-ICft) 

AN-0 0D0 

DO 322 II=*ICA,ICB 
X(II )=XC 1 )+N*ANXDELH*TAtFA/CBA 
Y(II)=Y( 1 )-N*DELH 
AN-AN+1.0DO 

E n*m+i 

AN-0 0D0 

ICD- ICB-ICA+1 

DO 323 J-JTC, JLC, JINT • 

X(J ) = X C 1 ) 

Y(J)»(1 0D0-AN/DOTH )*THICK 
DO 324 JJ-1. ICD 

X C J+ J J ) -X ( J )+ J JXEZERO.' C 2X C CBA+ 1 . 0D0 ) ) 

Y( J+JJ ) *Y ( J ) 

DO 325 K»1,KT 

XCJ+JJ+K J-XC J+JJ )+KXAPS/(2.0D0*KT) 

Y( J+JJ+K )«Y( J ) 
l AN-AN+1 0DO 
RETURN 
END 

SUBROUTINE F0RP1AF 
inPLICIT REAL&8 ( A-H - 0-2 ) 

DIMENSION B(3).C(3)/EA(3,3) 

COmON/PCRD/X ( 2O0 ) , Y ( 200 K EZERO , BUIH .THICK - TALFA . APS 
COnnON/NAFO/A C 200 , 50 ) - F ( 200 . 1 ) . H ( 2O0 ) . TS ( 200 ) , TSAT , CK . NEL . NSNP . h- 
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*H, I COR I 200.3) 

DO 101 N* 1 , NSNP 
F(N, 1 WO.ODO 
DO 102 Mh»1 , MEAN 
hCN, P1A )‘0. 0D0 
CONTI HUE 

DO 103 IEL = 1 , NEL 
IA-ICOR' IEL, 1 i 
IB* I COR*. IEL, 2) 

IC«ICOR( IEL , 3 ) 

B(1 IB J-VUC) 

E(2)=V(IC)-Y(IA) 

B<3 »-Vt Irt ;-V( IB ) 

C(1 )«X<IC>-X(IB > 

C(2J«XlIA)-X<IC > 

C(3)=X(IB)-XCIA) 

EL-D90RTCC(3 3 )**2 ) 

AS'DABS ( (B( 1 )*C<2)-B(2)*C(1 ) )/2 0D0 ) 
HC-HC IEL >/CK 
DO 104 J-l,3 
JJ-ICORCIEL, J) 

DO 105 K a 1 . 3 
KK-ICORC IEL , K ) 

EA< J,K J = C B C J)*8CK)+C( J)*C(K) )/(4XAS) 

IFCHC EG 0 OD0) GO TO 106 

HEL-HC*EL/6 ODO 

IF( J EQ 3) GO TO 106 

IF(K . EQ 3) GO TO 106 

IFCJ.EQ K) GO TO 107 

EA(J,K)=EA( J,K)+HSL 

GO TO 106 

EAC J,tC )«EA( J.K 3+2XHEL 
IF(KK.LT.JJ) GO TO 105 
NU-KK-JJ+1 
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h< JJ ■ Nt-J ‘ J A( JJ , NIJ )+EA( J • K ) 

0 CONTINUE 

c continue 

FE-HCtTSC IEL):*EL/2 ODO 
F‘. IA. 1 »*FCIA. 1 )+FE 
F ( ID, 1 '’F(IE-1 >+FE 
G CONTINUE 
RETURN 
END 

SUBROUTINE EANDEC ( NEQ . NAXB , NU£C ) 

IMPLICIT REALMS ( A-H , O-Z ' 

COMMON ' PCRD'-X < 200 ) , V t 200 ) . EZERO , BUIM , THICK . TALFA . APS 
COnnON/NAFO/A ( 200 , 50 ) . F ( 200 . 1) , H < 200 ) , TS ( 200 ) , TSAT, CK , NEL - NSNP , N~ 

* 

XN, ICOPf 200, 3 ) 

LOOP -NEQ- I 

DO 201 1=1 . LOOP 

riB«i+i 

riB«MIM0< I+flAXB-1 , NEQ ) 

DO 201 J-MB.NB 
L-J+S-MB 
D«ACI,L)/A(I, 1 ) 

DO 202 Pin-l.NUEC 

K! f(J.mm)-fcj,mm>-d#f(I.mm) 
mm-minocmaxb-l+i .neq-j+i ) 

DO 201 K-l.MM 
NN*L+K>1 

it. A( J,K )-A( J.K )-D*A(I , NN ) 

DO 203 I - 1 , NUEC 

F(N£Q, I )-F«NEQ, I )/A(NEQ. 1 ) 

DO 204 I =2, NEQ 
J-NEQ-I+1 

K-MINOC NEG-J+1 , P1AXB ) 

DO 204 MM -1. NUEC 
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DO 205 1 = 2, K 

ME' 0 J+L-l 

F < J • MM > -F ( J , MM >-h ( J , L < JtF ‘ MB , MM ) 
F f J , MM '*F( J,MM)/A(J, 1 ) 

RETURN 

EHD 
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